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ABSTRACT 



A method is developed in terms of the Cauer Second Form 
representation of continued fractions as a means of designing 
linear single-input single-output (SISO) control systems. 
Optimal closed loop solutions corresponding to a set of 
prescribed eigenvalues are obtained through minimization of 
a quadratic performance index. The partitioning method of 
the Cauer Second Form for system simplification is presented 
with a simplified inversion technique for the reduced order 
system. 
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I. INTRODUCTION 



The purpose of this research was to develop an algorithm 
for obtaining optimal closed loop solutions corresponding 
to a set of prescribed eigenvalues for single-input single- 
output (SISO) control systems. It was desired that the 
algorithm be adaptable to digital computer techniques and 
unrestricted by system order. 

The Cauer Second Form for system dynamics representation 
was chosen over other alternatives because of the regular 
pattern of the state and output matrices, and the method 
of linear system simplification. 

In Chapter II, several basic properties of both Cauer 
First and Second Forms are presented from the theory of 
continued fractions. A simple and efficient algorithm is 
also developed for inversion of the continued fraction in 
either form, independent of Routh’s algorithm. 

In Chapter III, the method of linear system order 
reduction based on the Cauer Second Form is amplified. 

The emphasis on this area was primarily to elucidate the 
various methods previously employed for. system simplifi- 
cation. 

The original theoretical work of this thesis is pre- 
sented in Chapter IV. The objective was to obtain closed 
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loop solutions corresponding to a prescribed set of eigen- 
values. While minimizing a certain cost function, which 
met desired system characteristics. It is shown, by examples, 
that the derived algorithm is equally capable of handling 
systems with multiple and/or complex, as well as, unique 
sets of real eigenvalues. 

The final chapter. Chapter V, presents a discussion of 
results and suggests areas for future study. 
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II. PROPERTIES OF CAUER FIRST AND SECOND FORMS 



A. CLOSED LOOP SYSTEM IN CAUER FIRST AND CAUER SECOND 
FORMS 



Consider the closed loop transfer function given by: 



Y(s) 

U(s) 



n-1 

Z 

L=0 



n 



s 



+ 



n-1 

Z 

i = 0 






a . s 
1 



( 2 - 1 ) 



with block diagram as given in Figure 2.1. Equation (2-1) 
can be expanded into the Cauer Forms of continued fractions 
as follows. 

1 . Cauer First Form 

a. Arrange the numerator and denominator poly- 
nomials in descending order. 

b. Perform continued division. 



Y(s) 

uTs7 



b,s +b„s + 
n-1 n-2 



* ^n-l * ^n-2 



. . . + b, s + b 
1 o 

s^ ^ + ... + a, s+a 
1 o 



( 2 - 2 ) 



hj^s + 



^2 



hsS + 



^4 " 



(2-3) 
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Figure 2.1. Block Diagram of an nth Order Linear System 
for Direct Programming (Bush Form) 




2. Cauer Second Form 



a. Invert the numerator and denominator and arrange 
the polynomials in ascending order. 



Y(s) 

uT?T 



a + a, s + 
o 1 

b + b^s + 
o 1 



n-2 . n-1 . n 

.a„s +a,s +s 
n-2 n-1 



■^^n-2® ^N-1® 



n-1 



(2-4) 



b. Perform continued division. 



Y(s) 

UT¥T 






^3 " 






(2-5) 



or 



1 




^3 ^ 



( 2 - 6 ) 



Block diagrams of both systems are shown in Figures 2.2 and 

2.3. 
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Figure 2.2. Block Diagram Representation of an Nth Order System 
(Cauer First Form) 




M 
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Figure 2.3. Block Diagram Representation of an Nth Order System 
(Cauer Second Form) 





B. PHYSICAL INTERPRETATION OF DOMINANT TERMS 
RESULTING FROM CONTINUED FRACTION EXPANSION 

It is known that the most dominant terms in equations 
(2-3) and (2-5) are the first quotients, h^s and h^, re- 
spectively. A meaningful interpretation for these terms 
can be accomplished by applying the initial value and final 
value theorems. Letting Y(s)/U(s) = F(s), by an asymptotic 
expansion approximation: 

1. For Cauer First Form 



lim f(t) ~ lim sF(s) ~ 



(2-7) 



t->o 



S-voo 



and 



lim f(t) ~ lim sF(s) ~ h 2 + h^ . 



( 2 - 8 ) 






S“>0 



2 . For Cauer Second Form 



lim f(t) ~ lim sF(s) ~ h 2 + h^^ 



(2-9) 



t-»-o 






lim f(t) ~ lim sF(s) ~ h^. 



( 2 - 10 ) 






s->o 



+ 



lim f(t) must exist. 
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Equations (2-7) and (2-10) are of considerable interest since 



they involve the dominant term, h^. The implication is that 
the Cauer First Form emphasizes the initial or transient 
response of the system; whereas, the Cauer Second Form em- 
phasizes the final or steady state response of the system. 

In general, the quotients lower in position ip the continued 
fraction expansion have less influence on the performance of 
the system as a whole, (h^ has less influence than h^, where 
i<j). Because many systems must meet a set of steady state 
conditions, the Cauer Second Form will be used for the 
prescribed eigenvalue problem. 



C. CONTINUED FRACTION INVERSION 

The theory of continued fractions was first associated 
with Routh ' s Algorithm by Wall in 1945, [1] and [2]. The 
following year Frank [3] extended and modified Wall's work 
to include complex coefficients. Both, however, applied 
Routh *s algorithm only to continued fraction expansions, 
not to the problem of inversion. 

In 1969, Chen and Shieh [4] developed an algorithm 
method for converting a continued fraction into a rational 
fraction of two polynomials. Their method, which makes use 
of Routh *s algorithm, is presented below. 

If the elements, h^, are known for any continued 
fraction, then the state and output equations can be 
written immediately from Figures 2.2 or 2.3. 
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"^l" 






^2 




^2^1 ^4^^l"'’^3^ ^6^^1^^3^ ■ ' '^2n^^l^^3 ^ 


^3 




^2^1 h^^(h^+h2) hg (h^+hg+hg ) . . . 

• * * « 






^2^1 ^6 ^^l'^^3''’^5^ ‘ ‘ ‘^2n^^l’^ ‘ ' ''’^2n-D 



1 

N 

M 

1 




1 


^2 




1 


CO 

N 


+ 


1 






• 






« 

• 


1 

c 

N 

J 




1 



z = Hz + Dr 



C 



*■^2 ^4 ^6 






z 

n 



( 2 - 11 ) 



( 2 - 12 ) 



( 2 - 13 ) 



C = L z 



( 2 - 14 ) 
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The coefficients of the characteristic polynomial, 

|sI-H1, become the first row elements of the required Routh 
array. The next sequence of steps in determining the (j+2)th 
row of the Routh array is to successively let 

h. = 0 
D 

and 

hj^.^ = 0 (2-15) 

for j s[l , 3 , 5 , . . . 2n_l] 

and evaluating the remaining ( n-k)x^n-k) matrix, where 
k = (j+l)/2, up to k=n-l, i.e., for n arbitrary and j=l; 
the 3rd row of the Routh array becomes (after h^ and h 2 are 
set equal to zero) the coefficients of 






where the (n-l)x(n-l) matrix is 



^ 4^3 ^ 6^3 



’^2n^3 



Va hg(h3+h3) h2^(h3 + h5) 



V3 



■^2n^^3'^‘ • ‘^2n-l^ 



(2-16) 



This process repeats until the system state matrix is 
reduced to a single element, yielding the (2n-l)th 

row in Routh ’s array. It is observed that each successive 
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odd numbered row contains one less element than it’s 
predecessor. By inserting leading zeros in the 3rd, 5th, 
..., (2n+l)th row, the matrix, P, is formed. 



3rd 


Pll 


P12 


P13 


1 


5th 


0 


P22 


P23 


1 


7th 


0 


0 


P33 


1 



(2-17) 



(2n+l)th 



0 



0 0 



1 



The matrix, P, is the linear transformation matrix required 
to obtain a linear system in Cauer Second Form from phase 
variable (canonical) form. Continuing, the second row of 
the Routh array is obtained from the output matrix, L, 
and the above transformation: 



c = Lz (continued fractions) 

z = Px (linear transformation) 

y = Cx (output equation, phase variable form) 

Therefore , 

C = L P . (2-18) 

C is an (Ixn) vector whose elements are the seoond row of 
the Routh array. 

Consider the Routh array as an (n+l)x(2n+l) matrix with 
typical element ♦ The quotients, h^, of the continued 
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fraction expansion can be expressed as: 



hi = r 



^il 



i+1,1 



(2-19) 



From this relationship and knowledge of how the Routh 
array is generated, the remaining even numbered rows of the 
array can be found. The transfer function as a ratio of two 
polynomials is written as: 



T(s) 



n 

Z 

i=l 

n+1 

Z 

L=1 





s 



i-1 



( 2 - 20 ) 



Chen and Shieh [4] contend that this method is the 
easiest in attaining the inversion. The author disagrees 
and presents a simpler iterative method based on the in- 
version technique for the Generalized Cauer Form given by 
Goldman [5]. The method is equally suited to both Cauer 
First and Cauer Second Forms, requiring no prior knowledge 
of Routh 's algorithm. Assuming all h^*s are known, and 
non-zero, in equation (2-3) or (2-5), let: 



a 



i 



h 



2i-l 



( 2 - 21 ) 



b 



i 




( 2 - 22 ) 



for i s[l,2,....,n]. 



20 



1 . Inversion of Cauer First Form 

Initialize two (n+lxl) vectors C and D. 



c = 


[Cq C 2 . . . 


n 


(2- 


23) 


D = 


[do d2 . . . 


... d ] 
n 


(2- 


24) 



to all zeros, except: 

c = b 
n n 

d T = a X c 
n-1 X n 

d = 1. 
n 



( 2-25) 
(2-26) 
( 2-27) 



The following set of equations are first solved for i=l. 



C • • 

n-i+] 



b . 

n-r 



X d 



n-i+j 



+ 



c 



n-i+j 



(2-28) 



d /•,-,\,«=a .xc .,.+d /•.-.N.., (2-29) 

n-(i+l)+] n-i n-i+3 n-(i+l)+]’ 

where je[0,l,2,...,i] are substituted in ascending order, and 
(2-28) is solved before (2-29) for each value of j. Now, 
let i=2 in equations (2-28) and (2-29) and repeat the same 
procedure. The index "i" is incremented until i=n-l, and 
(2-28) and (2-29) are solved as before over the appropriate 
range of the index "j". The final vectors, C and D, contain 
elements which are the coefficients of the numerator and 
denominator polynomials, respectively, of the transfer 
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function (or driving point impedance function): 



T(s) = 



n 

I C.s 
i = l ^ 



n -1 



n 

Z d.s 

j=0 3 



n-3 



Example : 



T(s) = 



10s^+171s+36Q 

s^+71s^+702s+720 



By continued fraction division: 



10s^ + 171s + 360 |s^ + 71s^ + 702s + 720^ 

s^+17 . ls^+36s 



,1s 



53.8s +666 s+ 720 



53.9s +666S+720 10s +171s 



+ 36o( 



10 



53.9 



10s^+123.562s+133.58 
47.438S+226 .42 



1 .1855 



47 .438S + 226 .42 |53 . 9s ^+66 6s + 7 2o( 

53.9s^+257. 258s 
408 .742S+720 



1.1362s 



408.742S + 720 /47 . 438s + 226 . 42\ .115 

47.438S+82.79 

143.63 



(2-30) 



(2-31) 



( 2-32) 
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143 .63 



'408.742s + 720C 2.8458s 

408.742s 

720 



720 [14 3 . 6 3(. .1995 

* 143.63 
0 



Therefore, the transfer function in the form of equation 
(2-3) is: 



T(s) = 



Is + 



.1855 + 



1.1362s + 



.115 + 



2.8458s + 1 



,1995 



(2-32) 



with 

hi = 

h2 = .1855 
hg = 1.1362 



h^^ = .115 

hg = 2.8458 

h^ = .1995. 
0 



(2-33) 



Now, using equations (2-21) and (2-22) j 



23 



.1 



.1855 



ai = 



bi -- h2 



1.1362 



= hg = 2.8458 



.115 

.1995 



(2-34) 



From equations (2-25) through (2-27) 



= .1995 

= ^ 3 X 0 ^ = 2.8458( .1995) = .568 

CI3 = 1. 



(2-35) 



Substituting i=l into equations (2-28) and (2-29) for j=0: 



= b 2 Xd 2 +C 2 = .115(.568) + 0 = .0653 



d^^ = a 2 XC 2 + d^ = 1. 1362( . 0653) + 0 = .0742 



for j = 1 : 



C 3 = b 2 Cd 3 +C 3 = .115(1) + .1995 = .3145 



d 2 = a 2 XC 3 +d 2 = 1 . 136 2( . 3145) + .568 = .9353 



At this point, j=i, therefore increment index "i" 
for i=2, j=0: 



= b^xd^^+c^^ = .1855(.0742) + 0 = .0138 



do = a^xc^+dp = .K. 0138) + 0 = . 00138 
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for j=l: 



C 2 = b^xd2+C2 



d 



1 



a^xC2+d^ 



.1855(.9353) + .0653 = .2388 
.K.2388) + .0742 = .0981 



for j=2: 

c'3 = b^xdg+Cg = .1855(1) + .3145 = .5000 

d2 = a.^xC^+d^ = .K. 5) + . 9353 = . 9853 

Now, at the point where j=i, and i=n-l, the transfer 
function is: 



T(s) = 



.0138s -<-. 2388s+.5 

.00138s^+.0981s^+.9853s+l 



(2-36) 



Multiplying numerator and denominator by 1/d^ = 720 yields: 



_ 10s^+171s+360 

s^+71s^+702s+720 



(2-37) 



2 . Inversion of Cauer Second Form 

Initialize the two (n+lxl) vectors C and D, (2-23) 
and (2-24), to all zeros, except: 



n 




(2-38) 


n-1 


= 1 


(2-39) 


“n 


a xc 
n n 


(2-40) 
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The following set of equations are first solved for i=l. 



n-l+j n=i n-i+j n-i+]+l 



(2-41) 



d /“inNi* — a • "a c ^d *it 

n-(i+l)+] n-i n-(i+l)+] n-i+1 



(2-42) 



where jE[0 ,1, 2 , . . . . ,i)^ are substituted in ascending order, 
and (2-41) is solved before (2-42) for each value of j. 



Next, find d^ according to: 



d 

n 



a 



n-i 



X c . 
n 



(2-43) 



Now, let i=2 in equations (2-41) and (2-42) and repeat the 
same recursive procedure. The index "i" is incremented by 
one until i = n-1, and for each value of i, (2-41) and (2-42) 
are iteratively solved over the appropriate range of the 
index "j". The resulting elements of C and D are the 
coefficients of the numerator and denominator polynomials 
of the transfer function: 

n-i 



T(s) = 



n 

E c.s' 
L=1 



n 

E d . s 

:=o " 



(2-30) 



n-1 



Using the sample example, (2-31); 



T(s) 



for j=i, 



10s^+171s+360 

s^+71s^+702s+720 

*^n-i+j+l = 0. 
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Place the numerator and denominator terms in ascending 
order, invert, and perform continued fraction division: 



360+171s+10s^ 



/720+702s+71s^ + s^ 
720+342s+20s^ 

360s+51s^+s^ 



2 



360s + 51s + s^ /360 + 171s + 10s^'- 1/s 

360+51s+s^ 

120s+9s^ 



120s + 9s 



2 



/360s+51s^ + s^^ 3 

360s+27s^ 

2 3 

24s ts'^ 



24s+s^ /l20+9sC 5/s (2-44) 

120+5S 
45 



4 /24 + s(. 6 

24 

s 



S /4\ 4/s 

4 
0 



The transfer function (2-31), in the form of equation (2-6) 
is : 



T(s) 



1 



2 + 



s 



1 + 



3+ 



27 



s 



5 + 



6 + 



s 



4 



(2-45) 



where 

hi = 2 h^^ = 5 

h2 = 1 h3 = 6 

hg = 3 hg = 4 (2-46) 

For the inversion process, using equations (2-21), (2-22) 
and ( 2-46 ) ; 

= 2 bi = h^ = 1 

a2 = hg = 3 = hj^ = 5 

ag = hg = 6 bg = hg = 4 (2-47) 



and from equations (2-38) through (2-40), 




dg = ag X dg = 6(4) = 24. ( 2-48) 

Now, substituting i = 1 into equation (2-41) and (2-42), for 
j = 0 : 



C 2 = b 2 + 63 = 5 + 4 = 9 

df = a 2 X Of + d 2 = 3(0) +1=1 (2-49) 
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for j=l; 



clg = b2 X dg + 0 = 5(24) = 120 

X ^2 + dg = 3(9) + 24 = 51 

and from equation (2-43): 

dg = a2 X Cg = 3(120) = 360 
for i=2, j=0: 

= C2 = 1 + 9 = 10 

(2-50) 

d^ = ag^ X Cq + dg^ = 2(0) + 1 = 1 

for j=l; 

62 = bg_ X d2 + Cg = 1(51) + 120 = 171 

di = ai X Cl + d2 = 2(10 + 51 = 71 

for j=2; 

dg = bg^ X dg + 0 = 1(360) + 0 = 360 

dg = ag_ X Cg + dg = 2(171) + 360 = 702 

and from equation (2-43): 

dg = a^ X Cg = 2(360) = 720. 
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Therefore , 



[0 


10 


171 


360] 




[1 


71 


702 


720] , 


(2-51) 



and the transfer function realized from equation (2-30) is: 



T(s) = 



10s +171S+360 
s^+71s^+702s+720 



(2-52) 



which is the same as (2-31). 

This completes the development of the continued fraction 
inversion algorithms from Cauer First and Second Forms. This 
iterative procedure is easily seen to be computationally 
much simpler than Chen and Shieh's method. First, it does 
not require the need to find the H matrix, (2-11); and 
second, it does not necessitate finding the coefficients 
of n characteristic polynomials of diminishing order. This 
method is solely based on equation (2-6), enumerating the 
inversion from bottom to top. As by-product, the entire 
Routh array appears in the intermediate steps as can be 
seen from the Cauer Second Form example: 



0 






720 

360 

360 

120 

24 

4 

1 



702 

171 

51 

9 

1 



71 

10 

1 



(2-53) 
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where rows 2(n-i)l and 2(n-i) are taken from the ith iteration, 
ilcCO ,1 , 2 , . . . ,n-l] . "i = 0" implies the rows come from the 
initialization of C and D. The last row, the (2n+l)th, is 
always the single element one. 

If a comparison is made between equations (2-20) and 
(2-30), it is observed that: 



= r 



2,n-i+l’ 



ie[l,2 , . . . ,n] 



d. 

3 



r^ j£[0,l,2,...,n3 , 



(2-54) 



where the c^*s and d ^ ' s are taken from the (n-l)th intera- 
tion under the index "i". 

It is also observed that if the quotients, h^'s, 
resulting from expansion into Cauer Second (First) Form 
are used in the inversion algorithm presented for Cauer 
First (Second) Form, then the c^ ' s and d^'s in the (n-l)th 
iteration represent the transfer function coefficients in 
reverse order. This is shown using the preceding example. 

From equation (2-46); 

- 2 h^ - 5 

hj = 1 hj = 6 

hj = 3 hg = 4 



and equation (2-47); 



ai - hi = 2 

a2 = hs = 3 
■^3 " ^5 ■ ® 



= h. 
= h. 



1 

5 

4 
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Now, using the inversion scheme for Cauer First Form, from 
(2-25) through (2-27); 



= ‘’S = ■* 

d, = = 6(4) = 24 (2-55) 



Making the substitution, i=l, in equations (2-28) and (2-29), 
for j = 0 ; 

= b2 X d2 + = 5(24) + 0 = 120 

d^ = a2 X C 2 + d^ = 3(120) + 0 = 360 



for j =1, ( j = i) ; 



= ‘’2 >^3 ^ '"’3 



dj = a2 X 03 4 dj 



for i=2, j=0; 



= bj^ X di 4 



<^0 = * '^0 



5(1) +4=9 
3(9) + 24 = 51 

* 

1(360) + 0 = 360 
2(360) + 0 = 720 



for j=l; 



C 2 = bi X d2 + C 2 



di = ai X C2 + d^_ 



1(51) + 120 = 171 
2(171 + 360 = 702 
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I 



for j=2, ( j=i) ; 



Cg = X dg + = 1(1) + 9 = 10 

d2 = X Cg + d2 = 2(10) + 51 = 71 



and dg remains 


unchanged , 


equal 


to 1 . 




Therefore ; 










•C = [0 


360 


171 


10] 




D = [720 


702 


71 


1], 


(2-56) 



and the transfer function should be; 

_ 360s^+171s+10 

T( s ) - 2 2 

720s'^ + 702s +71S+1 . (2-57) 

Since the h^'s from the Cauer Second Form were used in the 
inversion algorithm from the Cauer Frist Form, the vectors 
C and D require reversing non- zero elements, resulting in: 



C = 


0 

1 1 


10 


171 


360] 




D = 


[1 


71 


702 


720] , 


(2-51) 



and the correct transfer function is 

T(s) = (2-31) 

s +71s +702S+720 
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A digital computer program (FORTRAN IV) has been 
written for both Cauer First and Cauer Second Forms and is 
included as Appendix 3 with documentation. 
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III. LINEAR SYSTEM ORDER REDUCTION VIA PARTITION 
OF THE CAUER SECOND FORM 



A control system, in general, can consist of many tens 
or hundreds of elements. In such cases, the problems facing 
the engineer include: (1) too many variables to efficiently 

handle-, (2) the dimension of the system is too high to com- 
prehend; and (3) the modifications needed to meet required 
design characteristics are difficult to ascertain. A logi- 
cal approach is to seek procedural methods which reduce the 
order of the system to a manageable size yet maintain the 
basic characteristics of the full dimension model. 

A number of different methods for system simplification 
have been proposed for the reduction of high order dynamic 
systems to low order models of a more computationally or 
analytically tractable nature. The approaches used are 
quite different, but appear to fall into three main groups. 
The first is to ignore those modes of the original system 
which contribute little to the overall response. Davison 
[6] chose to neglect eigenvalues of the original system 
which are farthest from the origin, retaining only the do- 
minant eigenvalues and hence dominant time constants in the 
reduced model. The shortcoming of this technique is that 
many systems do not have any "dominant" roots [7]. 

Chidambara [8] essentially finds a reduced forcing function 
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so that the steady state values of the lower order model 
agree with those of the original system. The consequence of 
this method is that the approximate model give correct 
steady state values but incorrect time responses because 
the reduced forcing function does not excite the modes of 
the two systems in the same proportions [6]. Marshall [9] 
proposed the reduction of the state matrix by partitioning 
it and setting certain rate variables equal to zero in order 
to maintain the original steady state values. This techni- 
que, like Davison’s, is based on dominant roots and, there- 
fore, exhibits the same shortcomings. 

The second approach is to search in some manner for the 
coefficients of a set of differential equations of specified 
order, the response of which is sufficiently close to that 
of the original system when both are driven by the same 
inputs. Sinha and Pille [10] proposed a reduction technique 
based on the iterative application of the matrix pseudo 
inverse algorithm [11] to determine a model of specified 
order which minimizes the sum of the squares of the errors 
between the responses of the original system and the reduced 
order model to a given input. The main drawback of this 
method is that the objective function to be minimized is 
restricted to be the sum of the squares of the errors. Sinha 
and Bereznai [12] presented a method which minimizes a 
specified error criterion for a given reduced order model 
of the original system, based on the pattern- search algorithm 
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of Hooke and Jeeves [13]. Although this method provides 
more flexibility than that of Sinha and Pille, it generally 
requires considerably more computational time due to the 
poor convergence properties of the pattern-search algorithm. 

The third category involves application of the theory 
of continued fractions. Methods involving this approach 
have been developed by Chen and Shieh [14]. 

Sinha and DeBruin [15] and Fellows et al [16] have 
established the fact that among the methods previously 
mentioned, the approach by continued fraction expansion is 
generally the best for linear model simplification. 

A. SIMPLIFYING A TRANSFER FUNCTION 

The general nature of a control system is that of a low 
pass filter. Therefore, model simplification should con- 
centrate on the steady state aspects of the response with 
the transient portion given secondary consideration. As 
previously shown in Chapter II, the Cauer Second Form 
exactly characterizes these miens. 

Given the nth order original system transfer function: 
n-1 

S b^s^ 

T(s) = , (3-1) 

s^ + E a . s^ 
j=0 ^ 

where an mth order simplified model of the system (where 
m is strictly less than n) is desired, the polynomials in 
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equation (3-1) are rewritten in ascending order; 



T(s) = 






K ^n-1 

• • • • ]D 'O -} S 



S Q 3, ^ ^ * * * * 



.a ^+a , 

X-/ n-i 



(3-2) 



and expanded into a continued fraction: 



T(s) = 






h2 -H 



^3 ^ 



h2n-i ^ § • 

^^2n (3-3) 

An mth order simplified model is obtained by keeping the 
first 2m quotients of the expansion, omitting the remainder; 

T(s) = 1 



+ s 



'’2m-l * 



2m 



(3-4) 



and performing the inversion of the truncated fraction. The 
inversion technique presented in Chapter II can be used for 
this purpose . 



38 



Consider the seventh-order system, representing the 
control system of the pitch rate of a supersonic transport 
aircraft [10], described by its transfer function: 

_ 375000(s+. 08333) 

s"^ + 8 3. 6 3 5s® + 4097s^ + 7034 2s‘^ + 85 3 703s^ 

+ 2814271s^ + 3310875S+281250 (3-5) 

By continued fraction expansion: 

T(s) = 1 

9.00036 + s 

-.486286 + s 

-.036856+ s 



78496.2032+ ^ 

. 00071478 



(3-6) 



Suppose a second-order simplified model is desired. Equa- 
tion (3-6) has fourteen quotients. For the desired system, 
the first four quotients are kept with all others discarded. 
The truncated continued fraction becomes : 

T^Cs) = __1 

9.00036 + s 

-.486286 + s (3-7) 



-.036856 + s 

. 616185 
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and converted into transfer function form; 



.1299S+. 01105 
s^+1.14644s+. 09941 



(3-8) 



The block diagrams of equations (3-5) and (3-8) in the 
Cauer Second Form are shown in Figures 3.1 and 3.2, respec- 
tively. The unit step and impulse responses of the original 
and simplified systems are compared and shown in Figures 
3 . 3 and 3.4. 

B. STATE EQUATION SIMPLIFICATION 

The method of system simplification just presented is 
especially advantageous when converted into state space 
form. In Figure 2 . 3, a name for each state variable is 
given after each integrator, shown in Figure 3.5, from 
which the state equations and output equation can be 
directly written. 



— — 1 










^1 




^2^1 ^4hl ^6^1 ***^2n^l 




^1 


^2 




^2^1 ^6^^1^^3^ * ’ '^2n^^l"''^3 ^ 




^2 


^3 

• 




h2h^ hj^(h^+h3) hg(h^+h3+h5)...h2^(h^+h3+h5) 

• • • • 




^3 

• 

• 


• 

• 

_^n 




• • • • 

^2^1 \(V^3^ ^6^^l''^3'^^5^'- *^2n^^l''^3''"^2n-l^_ 




• 

z 
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9.0003 

Figure 3.1. Block Diagram of 7th Order System from Continued Fraction 
Expansion 



6162 
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Figure 3.2. Block Diagram of 2nd Order System from Truncated Continued 
Fraction Expansion 
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Figure 3.3. Unit Impulse Response Comparison of 
7th and 2nd Order Systems 



cro 
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Figure 3.4, Unit Step Response Comparison of 
7th and 2nd Order Systems 




3 
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Figure 3.5. Block Diagram of Ntli Order System Indicating 
State Variables 







1 

1 



1 



+ 



r 



and 



1 



(3-9) 



C - [h2 h0 • ♦ • *^2n^ 





(3-10) 



z = Hz + Dr , and C = Lz 



(3-11) 



Simplification of the equations in (3-11) can be achieved 
by partitioning of H, D and L, as indicated in Figure 3.6. 
The resulting mth order system becomes: 



Z" : 

~P 


’ ip 






where : 








II 


hjhi 


• 


^2m^l 

^2m^^l'’’^3^ 






h.Chi+hj) 


^2m^^l'^*** 



(3t12) 



(3-13) 
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Figure 3 #6, Nth Order System Simplification to Mth Order 
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D- 

~P 



1 

1 



1 



and = h z ‘ , 

p ~p ~p’ 



(3-14) 



(3-15) 



C ! 

P 



[h2 



"pi 

"P2 



z 

pm 



(3-16) 



As an example, consider the seventh order system described 
by the transfer function: 



C(s) 

R(s) 



1441.53s^+78319s^+525286 .125S+607693 .25 
s^+112.04s®+3755.92s^+39736.73s^ 



(3-17) 



+ 363650 . 56s^+759894.19s^+683656 .25s+ 617497. 375 



Arranging the polynomials in ascending order and expanding 
into the continued fraction yields : 
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Figure 3.7. 7th Order System Reduction to 2nd Order 
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Figure 3.8. Unit Step Response Comparison of 
7th and 2nd Order Systems 
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Figure 3.9. Unit Impulse Response Comparison of 
7th and 2nd Order Systems 




9rP 



C( s ) _ 

RTiT 



1 



s 



1.016133 + 

4.054112 + s_ 

-.067134 + s 



(3-18) 

595660.646 + s 

.0000757 

and equations (3-9) and (3-10) are formulated from the 
quotients. A simplified model of second order is desired, 
therefore, the state and output equations are partitioned 
as indicated in Figure 3.7. The simplified transfer 
function is: 



Cp( s) 
RP(s) 



. 250367S + 1. 035264 
s^+.509768s+l. 051966 



(3-19) 



Unit step and impulse responses of the original seventh and 
simplified second order systems are shown in Figures 3.8 and 
3.9 respectively. It is observed that the results of ex- 
pressing the seventh-order system by a second order model 
are satisfactory. 

It should be pointed out that a stable transfer func- 
tion may produce an unstable simplified function because the 
method of continued fraction expansion approximation does 
not necessarily guarantee a stable system. 
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IV. DESIGN OF OPTIMAL LINEAR CONTROL SYSTEMS 
WITH PRESCRIBED EIGENVALUES 



A. INTRODUCTION 

Consider the control of a plant with dynamics described 
by a set of first order, time-invariant linear differential 
equations of the form 

X = Ax + Bu, (4-1) 

where x is the nth-order state vector, A is the (nxn) 
plant matrix, u is the scaler control and B is the (nxl) 
input matrix. The output is defined as 

y = Cx , (4-2) 

where C is the (Ixn) output matrix. 

A linear feedback control law is assumed, and of the form 




There are mainly two separate approaches in the deter- 
mination of the feedback control matrix, G*, corresponding 
to the system under consideration; 1 - optimal control and 
2 - modal. 

In the optimal control approach a performance index is 

considered which is to be minimized in the design of a 

All states are available or an observer or Kalman filter 
is used to obtain the unknown states. 
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system. Assuming a performance index can be defined that 
represents most of the design requirements, "the solution 
to the optimal control problem can usually be obtained only 
by numerical methods that yield solutions to only a parti- 
cular problem" [17]. If solutions are sought to more than 
one numerical problem, simple performance indices must be 
defined, which often do not satisfy many of the design 
requirements. Therefore, the choice of a performance index 
must fall somewhere between a realistic criterion and one 
that is mathematically tractable. 

A quadratic performance index will be considered as a 
criterion for designing linear systems, of the form 

oo 

J - h f Cx^ Q X + Ru^]dt , (4-4) 

o 

where Q is a diagonal non-negative definite (nxn) matrix, 
and R is a positive scalar. 

In the modal approach, G* is chosen so that the closed 
loop system achieves the prescribed eigenvalues. Equations 
(1) and (3) together yield 

X = (A + BG*)x . 

If Q and R are given in the optimal control approach, 
then the eigenvalues of the closed loop system are uniquely 
determined, which may not realize the required performance 
characteristics or desired degree of stability for the 
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system. Using the modal control approach, a feedback control 
matrix can be found that will give the system the desired 
eigenvalues. This matrix is usually not unique, and it is 
not possible in a practical sense to find one that is 
"better" than it’s predecessors, since a performance measure 
is generally not known that corresponds to a given feedback 
control matrix. Therefore, it is necessary to find a method 
for determining the matrix, G* , that simultaneously satisfies 
the desired eigenvalues and minimizes a given performance 
index . 

In addressing this problem, Chang [18] and Tyler and 
Tuteric [19] have applied the root locus method to single- 
input, single-output and multivariable systems, respectively. 
The method lacks a rational computational procedure for 
determining the elements of the weighting matrix, Q, to meet 
a set of prescribed closed loop eigenvalues. Anderson and 
Moore [20] presented a restrictive method whereby a set of 
eigenvalues may be located to the left of a line parallel 
to the imaginary axis in the complex plane . Chen and Shieh 
[14] presented a method using sensitivity analysis. 

Solheim’s [17] method of a diagonalized (decoupled) system 
becomes complicated when the system contains either complex 
or multiple eigenvalues . Systems that cannot be diagonalized 
add further to the complication of the method. 

The method developed here takes advantage of the 
properties of the Cauer Second Form, is approached in a 
simplistic manner, and is easy to implement computationally. 
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V 




B. TRANSFORMATION TO PHASE VARIABLE FORM 

Consider an nth order linear system of the form 

e = Se + Tf , (4-5) 

with output 



d = We 



(4-6) 



Silverman [21], et al., have shown that if the system is 
controllable, then there exists a non-singular transformation 
matrix which takes an arbitrary state variable system to 
phase variable (canonical) form. (see Appendix A) 



X = Ax + Bu 



y = C X , 



where 



A = 



nl 



1 0 0 

0 1 0 0 

0 1 

• • • 

0 

0 0 1 

“n2 “nn 



(4-1) 

(4-2) 
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0 

0 



B 



0 



1 



C = [^^e. B. B] . 

~ j. Z o n 



If the system is not controllable, the phase variable 
(canonical) form may still be obtained from the system 
transfer function 
n 



Y(s) 

U(s) 



^ 8 s^^-l 
i=l ^i® 



n 

n , „ 

s + 1 



i = l 



a . s 
ni 



i-1 



(4-7) 



Once the system is in phase variable form, Chen and Shieh 
[22] have shown that the equivalent system in Cauer Second 
Form is easily written as 

z = Hz + Vu (4-8) 

y = C*z , (4-9) 

where the two forms are related by a linear nonsingular trans- 
formation matrix, P, 

z = Px . (4-10) 

The matrix P is an (nxn) upper triangular matrix. 



57 



The performance measure under consideration becomes 



T~ T 

J = ^ / [z Qz + u Ru]dt 



(4-11) 



where 



~ -IT -1 

Q = (P -^) Q P 



From optimal control theory, the Hamiltonian 



H = ^[z'^Qz + u'^Ru] + p"^[Hz + Vu] , 



(4-12) 



where P is the set of Lagrange Multipliers (also called the 
costate or adjoint vector) . For the Hamiltonian to be 
globally minimized, assuming no bounds on admissible states 
and control values, it is necessary that 3H/8u = 0 and 
3^H/3 u^>0 . 



3H/3U = Ru + V^p = 0 



(4-13) 



implies 



u* = -R~Vp , and 



3^H/3u^ = R>0 , 



(4-14) 

(4-15) 



since R was defined as a positive scalar. Included in the 
necessary conditions for optimality are 



z = Hz + Vu* 



3H/9z = -p = Qz + H p 



(4-16) 

(4-17) 
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Combining equations (4-14), (4-16) and (4-17) yields a set of 
2n linear differential equations forming the canonical sys- 
tem in Cauer Second Form. 





-1 T 

H VR 




z 


~ T 

-Q -H 




1 

1 



(4-18) 



It remains to be shown that this form is useful in obtaining 
optimal closed loop solutions that correspond to a set of 
prescribed eigenvalues. 



C. SIMILAR EIGENVALUES 



Consider the 2nth order cononical system in phase 
variable form: 



-1 T 

A -BR ^B 




1 

1 


T 

-Q -A^ 




1 

1 



(4-19) 



It is known that this system possesses n eigenvalues with 
negative real parts and n eignevalues with positive real 
parte, and that they are located symmetrically about the 
imaginary axis in the complex plane [17]. The eigenvalues 
of the optimal feedback system 

X = (A + BG*)x (4-20) 
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are identical to those with negative real parts in the 
canonical system. It is possible, therefore, to focus 
on the canonical system where the dependence eigenvalues 
on the weighting matrices Q and R can be studied without 
solving the matrix Riccati equation. The problem is in 
determining a weighting matrix, Q, such that the system 
attains the prescribed set of eigenvalues. 

The canonical system in Cauer Second Form can be ob- 
tained from the phase variable form using the linear trans- 
formation 

z = P X , (4-10) 



where P is a nonsingular matrix. We have 



X 



P 



1— -1 i — 1 

A -BR 




X 


T 

"9 




1 

1 




(4-19) 



(4-20) 



(4-21) 
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{•ti I j 2! i S J N . <>a > 



z 




PAP'^ 


-pbr'^b'^p'^ 




z 














A 

P 




-(P"^)*^ Q P'^ 


-(p"^)'^a'^p'^ 




1 

<P^J 
1 



(4-22) 





r -1 T 1 








H -VR -^V 




z 




~ T 




/N 




-Q -H 




P 



-1 T 

A -BR 



-Q 



-A' 



H 

"9 



-1 T 

-VR 

T 

-H 



5 



P 0 

0 



(4-24) 
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I 



where each sub-matrix of M, N, and P are known to be (nxn) 
square matrices. It is easily seen that P is non-singular, 
and that 





P 0 




P ^ 0 


p p"^ = 










0 (P"^)'^ 




T 

0 P 



Therefore , 



P P'^ 0 



^p-l)TpT 



0 

I 



(4-25) 



P M P"^ = N 



(4-26) 



shows the similarity of the M and N matrices. Two similar 
matrices have the following properties: 

1. Their determinants are the same. 



det M = det N 



(4-27) 



62 



2. Their traces are the same. 

Tr [M] = TrCN] (4-28) 

k k 

E m. . = E n. . (4-29) 

i=l j=l 

3 . Their characteristic equations are the same . 

det[XI-M] = det[XI-N] = 0 , (4-30) 

where X is an arbitrary variable. Since their characteristic 
equations are the same, the eigenvalues of M and N must be 
identical. It is now known that the Cauer Second Form and 
phase variable system matrices are similar in that they 
possess identical eigenvalues. 

D. DEVELOPMENT OF THE PRESCRIBED EIGENVALUE PROBLEM 

Initially given is a known linear system with dynamics 
described by either a set of first order differential equa- 
tions or its transfer function. It is desired to find the 
optimal feedback control, u* , such that the performance 
measure 



J = hi f [ X 



^11 



0 

^^22 




+ 
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u^Cl]u] } dt 



(4-31) 



is minimized, where the eigenvalues of the optimal system 
are specified as 



^ 1 ’ ^ 2 ’ ^ 3 ’ • 



1. Evaluation of the State (H) and Linear 

Transformation (P) Matrices 
' — ■ - ■ — -- ■- 

Assume the transfer function of the known system is 



given 



n 



T(s) = 



I & .s 
i = l ^ 



i -1 



s + Z . s 

i = l ^ 



(4-32) 



The system matrices in phase variable form are: 



A = 



0 


1 


0 


0 


0 




0 


0 


0 


1 


0 


0 




0 




• 


• 




• 


B = 


• 


-a. 


• 

• 

“OCo 


# 

• 




• 




• 

• 

1 


1 


2 


3 




n 







(4-33) 



C — 8/5 .... 8 ] 

~ 12 n 



(4-33) 



From equation (4-32), the Routh array is formed: 
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1 



“l “2 

^2 0 




d 



1 




0 



0 



IX] 

0 



(4-34) 



In matrix notation, the Routh array becomes j ^ j where 



ieCl, 2 , . . . . , 2(n+l) ] 

j e[l, 2 , . . . . , n+1] , 

and elements 

’^2(n-k) + 3, k ' ^ 

’^2(n-k)-4, k ■ ° 
k e [ 1 , 2 , . . . . , n+ 1 ] . 



In general, the (2L+l)th row of the Routh array is the Lth 
row of P . 
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p 



31 



0 

0 



32 


^^33 


3 ,n 


51 


r* ^ ^ . . . . . 


-p ^ 


52 


5 ,n-l 




r’„ , 
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7,n-2 



2n+l,2 

'2n+l,l 



(4-35) 



P being an (nxn) upper triangular matrix, will always have 
an inverse, P which can be quickly and efficiently deter- 

mined by digital computer. The H matrix formulation becomes: 



H = P A P 



-1 



(4-36) 





■^11 


hi2 


^13 


^in 




^11 


^22 


^23 


^2n 


= 




h^^ 


h,., 


h^ 






22 

• 


33 


3n 

• 






• 


• 


• 

h 






• 


• 


n-1 ,n 




h. . 


h„^ 


h.,^ 


h 




11 


22 


33 


n,n 


The elements 


of the 


H matrix 


can also be ' 


easily 


and directly 


from the 


first column 



(4-37) 
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Let 



= ? 



^i,l 

i+1,1 



(4-38) 



The h^’s correspond to the quotients of the continued frac- 
tion expansion of the transfer function in (32). 



T(s) = 






^2 " 



^3 



N . 



(4-39) 



The H matrix then becomes as shown in Figure 4.1. The 
regular pattern of the elements enable the H matrix to be 
obtained by inspection once the h^'s have been determined 
from either (4-34) and (4-38), or (4-39). 

The matrices V, and C*, are easily obtained: 



V = P B = 




(4-40) 
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Figure 4.1. State Matrix, H, in Cauer Second 



( 4 - 41 ) 



^ “ ^^2 ^4 ^6 * • * • 

B. Evaluation of Q 

Q = (P"^)'^ Q P"^ 





«11 


^12 


Qi3 ... 


^In 




^12 


^22 


Q23 


52 n 






— 






II 


^13 


^23 

• 


^33 

• 

# 


< 53 „ 




1 — I 

c 

• O' 

J 


• 

• 

«n2 


• 


• 

V 



The canonical system in Cauer Second Form (4-23) has now 
been obtained, with the numerical values of the elements 
of Q still to be found. This system will be compared to 
a non-optimized system with the prescribed eigenvalues also 
in Cauer Second Form. The desired system in phase variable 
form is: 



X = A*x + B*u 
y = E X 

Formulation icito Cauer Second Form yields 
X = H*z + V*u 
y = E*z 



(4-43) 

(4-44) 



(4-45) 

(4-46) 
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In a "nonoptimized" system (Q and R set equal to 0), the 
canonical system appears as: 




where 



H* 



0 




(4-48) 



possesses the n prescribed eigenvalues with negative real 
parts and n eigenvalues with positive real parts, symmetric 
about the imaginary axis in the complex plane. 

By equating the characteristic polynomials of N and N* , 
(4-23) and (4-47) respectively, it is now possible to 
determine the elements of the Q weighting matrix. 



det[sI-N] = det[sI-N*] 



det 



sI-H 

Q 



-1 T 
VR -^V 

T 

sI+H^ 



(4-49) 



det 



sI-H* 

0 



sI+(H*)'^ 



(4-50) 
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E. DETERMINATION OF THE OPTIMAL CONTROL LAW 

Finding the elements of the weighting matrix, Q, is 
obviously a non-trivial matter for all but the lowest order 
systems. The method developed for obtaining these values 
is based upon a succession of matrix building blocks, which 
are individually computationally simple. 

Starting with the matrix, H, 

H = [h^j] (4-51) 

define a new matrix, T, 

T = Ct^j] (4-52) 

where 



t . . 

ID 


= h . . - h . . 

ID DD 




i i D 


(4-53) 




n 








^ii 


= Z t . . 

j=i+l 


• 


i n 


(4-54) 



The matrix, T, therefore, is an (nxn) upper right triangular 
matrix, where the diagonal elements are equal to the sum 
of all other elements in the same row. The next "building 
block" matrix, G, is defined by: 



G = [g.j], (4-55) 

g.,=1.0 (4-56) 

^ 1,1 
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for ie[l, 2 , . . . , n] , 



§i,2 = ■^j-2,j-l 
for jg;[l, 2, n-1] , 




n-j + 2 
Z 

k=i+l 







(4-57) 



(4-58) 



for jeC3, 4,..., x] and iE[l, 2,..., n-j+1], where the 
index j is held fixed for each summation over the index i. 

The matrix, G, is an (nxn) upper left triangular matrix 
characterized by the first column being all ones. One more 
matrix needs to be defined at this point. Let the matrix, 

W, be defined as: 

W = (4-59) 

where 



w . 

- 1 . 






= (- 1 ) 



X g 






for i, j and k e[l, 2, ..., n] 



(4-60) 



The matrix, is therefore a tridimensional array with each 
"level" an upper left triangular matrix. Examples of the 
matrices, T, G and W are shown below. Although not evident 
at this point, the G and W matrices will be used heavily in 
obtaining the values of the elements of the Q matrix. 
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1 



I 



From the linear transformation 



Q = Q , 



(4-61) 



it is observed that once the element Q.. is known, the 
remaining elements in the same row (column) can then be 
obtained through a process of 



T = 



1 — 1 
1 — \ 
•p 


"^12 


• 

• 

• 

• 

00 
1 — \ 
p 


''^In 


0 


t22 


t23 


'‘^2n 




0 


‘'^33 


■^3n 






0 

• 


• 






• 

• 


• 

• 



nn 



(4-62) 



G = 



1 — 1 


1 — 1 
t—\ 

P 


gl3 ... 


®l,n-l 


1 — 1 


t22 


§23 ••• 


®2,n-l 


1— 1 


^33 


§33 


0 



1 

1 



n-l,n-l 
0 0 



®n-2,3 

0 



0 

0 



'l,n 



0 

0 



0 

0 



(4-63) 
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linear combinations of previously determined values. The 
Q matrix is thus found in the following manner: 



“ “ Q]_]_ “ ~ Qt <7 “ “ ^5 ~ “ ~ “ - - - _ _Q^^_ _ 



•12 '^13 



14 



■In 



*^22 " ■ *^23 " ■ *^24 ‘^ 2 n " 



^33 “ “ *^34 



•44 ’ 



■3n 



•4n 



Q - - 

^nn 



where for an nth order system (n>l) 



(4-65) 



•11 



2n „ 2n « 

n - n (h.)^ 

1=1 ^ i=l 

2n « 

n (h.) 

i=3 ^ 



2nH. « « 2 + 

n [(H.H, )^-(h.h, ) ] 

1=3 >'i ^ ^ ^ ^ 



(4-66) 



* for n=l, = (H^H^)^ - (h 2 h^)^. 
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Q. . 

13 



j-1 P 



- Z 



k,n-l 



k=l n-l,n-l 



<5ik 



(4-67) 



for ii^j or j iiin , 



‘in 



n-1 

Z Q. . 

j=i 



(4-67) 



and 



^nn 



n 

Z 



n ^ 



n 



n 



i=l j=l 



Z Q. . + ( Z H. . ) 



•13 



i=l 



11 



( Z h. . ) 
j=l 



n n n n 

+ 2C Z Z - Z Z h. . (h. .-h. . )] (4-69) 

i=l j=l L=1 j=l 



The values of the diagonal elements of Q (other than Q]_q^) 
generally involve varying linear combinations of already 
determined values, more easily expressed in terms of the 
and W matrices, rather than the P and H matrices. To aid 
in determining the values along the diagonal the following 
labels are provided for G and W. 



"P. . 

i>3 



P 
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(4-71) 
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These two arrays will provide the coefficients of each "Q..x 

]c 

s " necessary. Assume, for example, that a 5th order 
system is being considered in equation (4-50). The results 
of the expansion of both sides of the equation results in: 



2i “ 2i 

a^ + Z a.s = b + Z b.s (4-72) 

° i=l ^ ° L=1 ^ 

where a^^ = 1 , b^ = 1 . By equating and b^ ^ for 

i£[l,2, ..., n] , it is possible to obtain Q... 



i.e. a = b (4-73) 

o o 

a and b are the coefficients of s'^. From (4-70), it is 
o o 

observed that the only element in the s*^ column that is 
non-zero is which appears in row 1. The 1 "indicates" 

that it is necessary to only look in "level" 1 of W, in the 
column corresponding to s*^. This yields only the single 
element Therefore, the coefficient of will be 

Sl5 == ”l51 • 



What remains to be determined are the coefficients of s^^ 
not involving the Because of the symmetry of the 

eigenvalues, of both the system being "optimized" and the 
system with the prescribed eigenvalues, the remaining 
coefficients (those not involving a Qj^j) 3.re easily 
determined from: 
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det[sI-H] X detCsI h'^] 



(4-75) 



and 

det[sI-H*] X det[sI+(H*)'^] , 



(4-76) 



which give: 



( Z a. s^)( Z (-l)^a. s^) (4-77) 

i=0 ^ i = 0 

and 



i ^ k i 

( Z a.*s )( Z (-1) a.*s ) (4-78) 

i=0 ^ i=0 ^ 



where k=i+l for n even and k=i for n odd. The indicated 
multiplication in (4-77) and (4-78) result in: 



n 

Z 

i = 0 



(-l)^a^s^^ 



(4-79) 



and 



n 

Z 

i = 0 



( 1 A^2i 

(.-1) a^'*s 



(4-80) 



respectively, where the same conditions are imposed on "k" . 

Returning to the 5th order example, the coefficients of 
s^ are -Oq and -cJq*. To obtain is now a matter of 

solving the equation: 



(g 



15 



W ) 
'"l51^ 



•11 



- = - 0 , 



(4-81) 
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(4-82) 



^0 - 



'll 



(s 



15 



g 



151' 



With known, it is a simple matter to obtain the 

remaining elements in the first row (and column) using 
(4-67) and (4-68). 

To find Q 22 ’ 3,11 the coefficients of s . s results 

in three separate ways: 



( 1 ) 



0 2 



S X s 



(2) s^ X s^ 

(3) s^ X s*^ 



(4-83) 



The multiplicand indicates which columns of G, (4-70), is 
of interest. Any non-zero element in G tells which level 
of W is of interest. The multiplier is the indicator for 
the column of interest in the array, W. Therefore, for 



§15 ^ ^^131 ^11 ^231 ^12 ^331 ^ 13 ^ 



for s^ X s^: 



(4-84) 
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‘5ll 
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Qi2> * 




2,4 ^ 


(^142 


'^21 


+ ^242 


« 22 > 


( 4 - 85 ) 
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^20 
for s X s : 

+ 

+ 

(4-86) 

Obtaining Q ^2 is now a matter of solving the equation given 
by the combination of (4-79), (4-80), (-84), (4-85) and 
(4-86) : 



Si, 3 ^ ^'^ISI ^11^ 



^ 2,3 ^ ^'^152 ^ 21 ^ 



S3 3 ^ 



153 ^31' 



Sis ^i31 Qj^) 



®j,3 ^15j ^ji^ “ ^1 " “^1* ’ 



(4-87) 



for Q 22 * 

Once Q 22 has been found, the remainder of the 2nd row 
(and column) can be obtained using (4-67) and (4-68). 
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respective 
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(4-88) 
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( 



I 



1 



I 



3 5 - 4 4 ^ 

( Z g. 2 W. . Q..) + ( Z g. Z W . Q..) 

j=l i=l i=l i=l 



* ‘ji ®j.l j/i3j " “3 = ■ 



(4-89) 






(4-90) 



The underscore in (4-88), (4-89) and (4-90) indicates the 
term that contains *^44’’ *^55’ respectively. 

The entire Q matrix is now known. By the inverse 

A,^ 

transformation of (4-42), 



Q = p'^ Q P . (4-91) 

Q along with A, B and C of (4-1) and (4-2) are used to 
obtain the matrix Riccati equation steady state solution. 

0 = KA + A^K - KBR"^b'^K + Q • (4-9 2) 



Once K has been determined, the optimal control law is 
given by 

u* = -R“^BKx • .^ (4-93) 

-%ir ^ ^ ^ 

Once u* has been determined, an inverse non-singular trans- 
formation can be performed to take the phase variable form 

* It is now known that G=-R ^b'^K in (4-20). 
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back into state variable form. (See Appendix A) 



F. ILLUSTRATE EXAMPLES 

1 . Odd Order System (Third) 
Given the system 



^1 




0 


1 


0 




^1 




0 


^2 


= 


0 


0 


1 




^2 


+ 


0 


_^3_ 




-13 


-19 


-7 




_^3_ 




1 



y = [1 0 0] 



X, 



X, 



(4-95) 



find the optimal control law, u* , such that the quadratic 
performance index (4-31) is minimized, where the eigenvalues 
of the system are specified as 



s 



1 



-3, 




-4, 
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I 

I 



1 

I 

I 

I 

J 



Forming the Routh array yields 



13 

1 

19 

-7/19 

30/7 

1/30 

1 

0 



19 

0 

7 

-1/19 

1 

0 



7 

0 

1 

0 



1 

0 



from which the third, fifth and seventh row are extracted to 
form P. 



P = 



19 7 1 

0 30/7 1 

0 0 1 



(4-97) 



and 



,-l 



1/19 

0 

0 



H=PAP 



-1 



-13/19 

-13/19 

-13/19 



-49/570 


1/30 








7/30 


-7/30 






(4-98) 


0 


1 








1 ) , H and V 


are calculated 




637/570 


-13/30 






-63/19 


9/7 




(4-99) 


-63/19 


-3 
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I 



V = PB = 



1 

1 

1 



( 4 - 100 ) 



The set of prescribed eigenvalues yields the system: 



0 

0 



1 

0 



0 

1 



-72 -54 -13 









- 0 ~ 




1 








^2 


+ 


0 








1 




__ 3 _ 







u 



( 4 - 101 ) 



and 



y = [1 0 0] 



72 

1 

54 

- 13/54 

115/13 

1/115 

1 

0 



54 

0 

13 

- 1/54 

1 

0 



X, 



From the Routh array 



13 

0 

1 

0 



1 

0 



( 4 - 102 ) 



( 4 - 103 ) 



85 



} 

I 

I 



I 



extracting rows 3, 5 and 7 yields the P* matrix: 



54 

0 

0 



13 

115/13 
0 



1 

1 






1/54 -169/6210 1/115 

0 13/115 -13/115 



H* =P*A* ( P* ) 



-1 





1 _ 


_ 0 0 


1 


are 


calculated , 




(4-104) 


-4/3 


676/345 


-72/115 




-4/3 


-2574/621 


1980/1495 




-4/3 


-2574/621 


-11254/1495 


(4-105) 



V* = P*B’ 



1 

1 

1 



(4-106) 



Formulation of equation (4-50) yields 





si - H 


< 

1 

< 

1 




sI-H* 


0 


det 


Q 


T 

sI+H^ 


= det 


0 


T 

Sl+H 
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det 



"s+13/19 


-637/570 


13/30 


13/19 


s+63/19 


-9/7 


13/19 


63/19 


s+ 3 


\l 


Qi2 


?13 




Q 22 


^23 


^13 


^23 


^33 



1 


1 


1 


1 


1 


1 


1 


1 


1 


s-13/19 


-13/19 


-13/19 


637/570 


s-63/19 


-63/19 


-13/30 


9/7 


s-3 



■3 + 4/3 


-676/345 


72/115 


0 


4/3 


s+2574/621 


-1980/1495 


0 


4/3 


2574/621 


s+11245/1495 


0 



0 0 
0 0 
0 0 



=det 



0 



s-4/3 -4/3 -4/3 



0 0 0 
0 0 0 



676/345 s-2574/ -2574/ 
621 621 

-72/115 1980/ s-11245/ 
1495 1495 



(4-107) 



It is desired to determine the values of Q. Brute force 
enumeration of the determinants and equating coefficients 
of like powers of s would eventually lead to the solution. 
Using instead the method developed, from H(not H*) evolve 
the T, G, and W matrices. 



T 



7;0000 4.4333 2.5667 
0 4.2857 4.2857 
0 0 0 



(4-108) 
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I 



1 7.0000 

1 4.2857 

1 0 



19.0000 

0 

0 



(4-109) 



and W = 



-1 7.0000 
-1 4.2857 
-1 0 



-19.0000 

0 

0 



(4-110) 



Now add the appropriate labels as in (4-70) and (4-71) 



for G: 



1 

2 

3 



s 

1 

1 

1 



s 

7.0 



0 

s 

19.0 



4.2857 0 



(4-111) 



for W: 



911 

912 
9is 



-1 

-1 

-1 



7.0 -19.0 

4.2857 0 



0 



0 



(4-112) 



Still necessary are the coefficients in (4-79) and (4-80): 



( E a-s^)( z (-l)^a.s^) = Z (-l)^a.s^^ 
i = 0 i = 0 i = 0 



(4-113) 
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I 



{ 



( 



I 

I 



f 



for n odd (n=3), k=i+l, 



( Z a.s^)( Z s^) = Z 

i=0 ^ i=0 ^ i=0 ^ 



(4-114) 



(s^+7s^+19s+13)(s^-7s^+19s-13) 

6 4 2 

= s - 11s + 179s - 169, 

and 



3 *3 • . , 

( Z a.*s^)( Z (-1)^ -^a.^s^) 
i=0 ^ i=0 ^ 



3 

Z (-1) 
i = 0 




(4-115) 



(4-116) 



(s^+13s^+54s+72)(s^-13s^+54s-72) 
= s^ - 61s^ + 1044s^ - 5184 . 



(4-117) 



Now, by equating coefficients of s^, it is possible to 
obtain From G, the only non-zero element in the "s*^ 

column” is 19, which corresponds to row 1, and therefore 
level 1 of W. From level 1 of W, the only non-zero element 
in the ”s^ column” is -19, with row label Q 2 _]_- The 
coefficient of as obtained from (4-81) is 19(-19) = 

-361. The solution for is obtained from equations (4-81) 

and (4-82); 
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-361 - 169 = -5184 

Qll = (-5184 + 169) / -361 
= 13.892 . 



From (4-65), (4-67) and (4-68), Q ^2 success 

obtained: 





(13.892) = -22.69 



^13 “ “^11 ■ ^12 ■ ■13-892 - (-22.69) = 8.798. 



The next power of s which results from expansion of (4- 

2 2 
is s . Equating coefficients of s , it is now possible 

~ 2 0 2 1 1 
obtain Q 22 * ® results from the products s xs , s xs , 

2 0 

and s xs . Therefore, from the development starting at 
(4-83), the equation to be solved for is: 



1 

( 2 g 

:=i 



j,3 



Z w.^. Q..) t ( 



i = l 



2 

I i 
j=l 



2 

i 2 ^ ^ 

3 i=l 



i 2 j ’’ji’ 



3 1 

( 2 g. i 2 

j=l 3’^ i=i 



W... Q. .) 

i3] ^31 



+ a. 




(4-118) 

ively 

(4-119) 

(4-120) 

107) 

to 



(4-121) 
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I 

I 

1 







Substituting known values, equation (4-121) becomes; 



®13^^11l'^ll '*■ ^211 ^12 ^311 ^13^ 

gl2(Wi2i + ^221 ^12^ ^22^^122 ^21 ^222 ^22^ 

+ §11^^131 ^11^ ^21^^132 ^21^ ^31^^133 ^31^ 

+ (4-122) 

19 .0C-K13 .892 ) - K-22.69) - 1( 8.798)] 

+ 7.0 [7.0(13.892 ) + 4 . 2 8 5 7 ( -2 2 . 6 9 ) ] 

+ 4.2857 [7.0(-22.69) + 4.2857 Q22^^ + l[ -19 ( 13 . 89 2 ) ] 

+ 1 [-19(-22 .69 )] + 1[-19(8.798)] + 179 = 1044 

Q 22 = 84.155 (4-123) 

From (4-68) 

*^23 ■ "'^21 " *^22 " “^12 " ^^22 

- (-22.69) - (84.155) = -61.465 (4-124) 
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is the next power of s obtained in the expansion of (4-107) . 

Equating coefficients of s , it is now possible to obtain the 

4 

equation from which can be found. s results from multi- 

2 2 

plying s xs . Therefore, from G and W: 
gllC-KQll) - - ^^^13^^ 

+ g2it-KQ2i) - 1(Q22^ " 

" g3it-l(Q3l) - 1(^32^ - 1^^33> 



02 = - 02 » , 



(4-125) 



where g-j _]_5 g 2 q S 3 I equal one. Therefore, (4-125) 

becomes 



3 

E 



3 

E 



Q . . - o^ 
i=l j=l 2 



-02« 



(4-126) 



Solving (4-126) for Q 33 ' 



■33 



•02” '^2 



3 

E 



Since 



Q* 

^in 



n -1 ~ 

^ Qik 

k=l 



i=l j=l 
i+ j 6 



3 - 

E Q. . 



(4-127) 



(4-68) 
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(4-128) 



n ^ 



Z 

k=l 




0 , 



and (4-127) simplifies to 




2 

Z 

k=l 







- 61-11-8.798+61.465 = 102.667 



Q is now entirely known. 



Q = 



13.892 

-22.69 

8.798 



-22.69 8.798 

84.155 -61.465 

-61.465 102.667 



From ( 4-91 ) : 



Q = 




5015 

0 

0 



0 0 

865 0 

0 50 



(4-129) 

(4-130) 



(4-131) 



(4-132) 



To find the optimal feedback control law, u* , it is necessary 
to solve the matrix Riccati equation (4-92), where Q, A, B 
and C are as given in (4-131), (4-94) and (4-95). The 
solution yields: 
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-[59 35 



6 ] 



(4-133) 



U' 



or u* = -5 9 - 35 X 2 - 6 

With G* as given in (4-133), x = (A + BG’*)x 
becomes 




(4-134) 



Comparing (4-134) with (4-101), it is observed that (A+BG*)= 
A*; therefore, the desired system has been realized with the 
set of prescribed eigenvalues. 

2 . Even Order Example (Fourth) 

Consider the linear system represented by the transfer 
function: 



T^(s) 



s +9s + 34 

s^+12s^+48s^+80s+48 



(4.135) 



the system eigenvalues are -2, -2, -2, and -6. 
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This system could be expressed in phase variable form (4-31), 
thereby obtaining the transformation matrix, P; the system 

'V 

can be realized in Cauer Second Form (4-36). Instead, in 
this example, it is desired to obtain the quotients which 
result from partial fraction expansion of the transfer func- 
tion. In Cauer Second Form, T(s) becomes: 



T(s) 




49 + 80s + 48s^ + 12s^ + s‘^ 



1 




s 




s 




s 




s 




s 




s 




(4-136) 



s 



h 



8 
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where the quotients, , 



are : 



= 1.41176 
= .505245 
h3 = -4.6287 
= -.627918 



hg = 23.07 

h_ = .146914 
6 

= -281.808 
hg = -.024241 



Substituting these values into the matrix in Figure 4.1 
yields 



.71328 


-.88647 


. 20741 


-.03422 


. 71328 


2.01997 


-.47261 


.07798 


.71328 


2.01997 


2.91670 


-.48125 


.71328 


2.01997 


2.91670 


6.35004 



Again, V is (and in all cases considered will be) a vector 
of all ones , 



V = 




(4-138) 



It is desired, as before, to find the optimal control law, 
u* , such that the quadratic performance index (4-31) is 
minimized, and that the optimal system realizes a set of 
prescribed eigenvalues. The eigenvalues are given as: 
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J 



I 

I 



-5 




- 2 , 



s 



2 



S 3 , s^^ = -6 + j3 



(4-139) 



With the zeroes of the transfer function (4-136) the same, 
the transfer function of the desired system with the pre- 
scribed eigenvalues becomes: 



T^^*(s) = 



s +9S+34 

s‘^ + 19s^ + 138s^ + 435s + 4 50 



(4-140) 



By continued fraction expansion of T*(s) the quotients 
obtained are: 



= 13.2353 
H 2 = .107635 
H 3 = -71.3206 
= -.088175 

from which H* is: 



= -1077.43 

= .004834 
b 

= 396.926 
Hg = -.024294 , 



H 



5 % 



1.42458 

1.42458 

1.42458 

1.4258 



-1.16702 

5.12168 

5.12168 

5.12168 



.06399 

-.28082 

-5.48975 

-5.48975 



-.32154 

1.41114 

27.58655 

17.94349 



(4-141) 
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Formulation of equation (4-50) yields in matrix notation: 



det = 


sI-H 


1 

T 




det 


sI-H* 


0 

T 




Q 


sI+H^ 






0 


sI+(H*) 



which is to be solved for Q. 

As before, we formulate the T, G and W matrices 
according to equations (4-53) and (4-54); (4-56), (4-57) 
and (4-58); and (4-60), respectively. 

This yields: 



T 



12.00 

0 

0 

0 



2.9064 

9.6614 

0 

0 



2.7093 

3.3893 

6.8313 

0 



6.3500 

6.2721 

6.8313 

0 



(4-143) 



1 


12.0000 


46.5882 


67.2941 


1 


9.6614 


23.1534 


0 


1 


6 .8313 


0 


0 


1 


0 


0 


0 



(4-144) 



the notation 1 indicates an (nxn) matrix with all elements 
unity . 
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i CD 



w 



•r 



-1 12.0000 
-1 9.6614 

-1 6.8313 

-1 0 



-46 .5882 
-23.1534 
0 
0 



67 .2941 
0 
0 
0 



(4-145) 



With the proper labels, (4-70) and (4-71) 
for G : 

2 1 



1 

2 

3 

4 



s 

1 

1 

1 

1 



12.0000 
9.6613 
6 . 8313 
0 



46.5882 

23.1534 

0 

0 



0 

s 

67.2941 

0 

0 

0 



(4-146) 



for W; 





.•r 


1 

s 


0 

s 







/ 




✓ 


Qil 


^-1 12.0000 


46 . 5882 


67 . 2941“ 






-1 9.6613 


23.1534 


0 




<5i3 


-1 6.8313 


0 


0 


-7 


«i4 


O 
1 — 1 


0 


0 


/ 



(4-147) 



2k 

and W contain all of the coefficients of Qj^jS resulting 
from the expansion of (4-142). Still to be found are the 
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coefficients of s not involving any Qj_j > the and 

Therefore, from (4-79) and (4-80): 

(s^+12s^+48s^+80s+48) (s^-12s^+48s^-80s+48) 

= s® - 48s® - 480s‘^ - 1792s^ + 2304, (4-148) 



Oq = 2304 
= -1792 
= 480 

and 



= -48 



a 



4 



1 



(s^+19s®+138s^+435s+450) x 
(s‘^-19s^ + 13 8s^-43 5s + 450) 



= s®-85s® + 3414s*^ - 65025s^ + 202500 , 



Qq* = 202500 
0 ^' = -65025 
02* = 3414 



03>^ = -85 

^ 4 * = 1 



(4-149) 



(4-150) 



(4-151) 
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r 



Determination of the Qj_j's results from equating 

coefficients of like powers of s from the expansion of the 

determinants in equation (4-142). The Q^^’s are obtained 

in a successive method (4-65) by starting with the s*^ 

coefficients, then continuing by equating coefficients of 
2k 

s , ke[l, 2, ... n-1], in scandent fashion. 

For s*^: 

from G, the only non- zero element is corresponding to 

row 1, therefore, level 1 in W. The only non-zero element in 
level 1 of W in the s^ column is Equating coefficients 

yields : 



^11 * = "O" 

Solving: 

; _ 202500-2304 

^11 " (67.2941)(67.2941) “ 

and from (4-67) and (4-68), 



'^12 


= Q2I 




-88.95327 


^13 


■ ^31 




48.14819 


«14 


= Q4I 




-3.40294 



(4-152) 



44.20803 (4-153) 



(4-154) 
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2 

for coefficients of s ; 



1 

Z 



j=l 




3 

Z W.». 
i=i ^ 2 : 




2 

z 

j=l 




z w.„. 

i33 



i = l 




3 

+ 2 g 

j = l 



j2 



J/i4j 



+a^* 



Solving : 



(4-155) 



Q 22 = 296.94152 (4-156) 

and from (4-67) and (4-68), 

«23 = «32 ■- - 263-70162 

Q 24 = Q42 = 55.71337 (4-157) 



for coefficients of 



4 

s : 



2 

Z 



j=l 




4 

Z 

i=l 



%lj 






3 

Z 

i=l 



W 



12j 




■f 



4 

2 g 

j=l 






2 

z w 

L=1 



i3j 





(4-158) 
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Solving : 



Q 23 = 351.24159 



(4-159) 



and from (4-6 8 ) , 



= Q ^3 = -135.68816 



for coefficients of 



6 



s 



4 

5 



j=l 




4 

Z W 
i=i 



ilj 





(4-160) 



(4-161) 



4 4 

Q = - E Z Q.. + 03 * - a (4-162) 

j=i i=i 13 ^ ^ 

i+ j /8 



-Q 



44 



- 2 Q . 

i=l 



+ a. 



- a. 



(4-163) 



= 120.37773 



(4-164) 



Q is now entirely known. 



103 







\ 



Q 



44.20803 

-88,95327 

48.14819 

-3.40294 



-88.95327 

296.94152 

-263.70162 

55.71337 



48.14819 

-263.70162 

351.24159 

-135.68816 



From equation (4-91): 



-3.40294 

55.71337 

-135.68816 

120.37773 

(4 





“ 200196 


0 


0 


0 


T'' 

? 9? = 


0 


63233 


0 


0 




0 


0 


2934 


0 




0 


0 


0 


37 



(4 



Substituting matrices Q, R, A, B and C into the Riccati 
equation, and solving, yields the optimal feedback control 
law , u* ; 



u* = 



-[402 355 90 7] 



X. 



X, 



(4- 



where G* = -[402 355 90 7]. 



(4- 



The optimal closed loop system, x = (A+BG*)x, 
in phase variable form is: 



165 ) 



166) 



167) 



168) 
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1 



I 



1 



\ 



I 

I 



I 







“■ 0 


1 


0 


0 




^1 


X2 




0 


0 


1 


0 




^2 


^3 




0 


0 


0 


1 




^3 


_^4_ 




-450 


-435 


-138 


-19 




^4 






— 













(4-169) 



which realizes the desired set of eigenvalues. 

3 . Higher Order Example (Seventh) 

The previous examples represent an odd and an even 
ordered system, illustrating the minor differences in com- 
putational procedures. It is observed that for low order 
systems, the calculations can be done, relatively easily by 
hand. Higher order systems require, laborious and tedious 
computations . Appendix C provides a digital computer program 
which yields the weighting matrix , Q , with the only required 
input being the transfer functions of the known and desired 
eigenvalue systems. 

The following seventh order example utilizes the 
results from the program given in Appendix C. 

Consider the linear system given by its transfer 
function : 



T^(s) = 



249.435788 

s'^ + 9 . Os^ + 40 .4s ^ + 116 .8s^ + 233.6s^ + 323.2s^ 



(4-170) 



+ 288.0s + 128.0 
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It is known that the characteristics of the desired system 
are such that its resulting transfer function is given by; 



T^*(s) = — 



249.435788 



s'+15,4s®+101.64s^+372.68s^+819 ,896s^ 



(4-171) 



+1082.26272S +793.659328S+249.4357888 



The diagonal Q matrix was determined from the computer pro- 
gram yielding: 



• 5 ll 


= 


45834.21273428 


Q22 


= 


89780.21841734 


^33 


= 


55970.39786199 


Q44 


= 


19170.7157376 


Q55 


= 


3959.33664 


^66 


= 


494.9776 


Q 77 


= 


33.68 



Q, R and the state and output matrices representing the 
transfer function in equation (4-170) in phase variable form 
were substituted into the matrix Riccati equation. The opti- 
mal feedback control law, u*, resulting from solution of the 
Riccati equation is: 
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u* = -121.435789 - 505.659319 

-759.06269 - 586.29596 x^ - 255.879974 

-61.2399914 Xg - 6.39999944 x^ , (4-173) 

where the ma-ferix G* is: 

G* = -Cl21. 435789 505.649319 759.06269 586.29596 

255.879974 61.2399914 6.39999944]. (4-174) 

The optimal closed loop system, 

X = (A + BG*) X 

y = Cx , 

expressed as a transfer function is: 



Y(s) 



249.435788 

s^+15 . 39999974s®+101. 6399914s^ + 



372.679974s^+819.89596s^+ 

1082. 26289s^+793. 659319S+249.435789, 



which realizes the desired transfer function with the pre- 
scribed eigenvalues. 

Figures 4.2 through 4.13 show the impulse and step res- 
ponses of the three previous examples, both before and 
after compensation. 
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a 
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Figure 4.2. Unit Impulse Response of the Uncompensated (1) 

and Prescribed Eigenvalue (2) Third Order Systems 
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Figure 4.3. Unit Step Response of the Uncompensated (1) and 
Prescribed Eigenvalue (2) Third Order Systems 
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Figure 4.4. Unit Impulse Response of the Compensated (1) and 
Prescribed Eigenvalue (2) Third Order Systems 



1 
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Figure 4.5. Unit Step Response of the Compensated (1) and 
Prescribed Eigenvalue (2) Third Order Systems 
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Figure 4.6. Unit Impulse Response of the Uncompensated (1) and 
Prescribed Eigenvalue (2) Fourth Order Systems 
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Figure 4.7. Unit Step Response of the Uncompensated (1) an 
Prescribed Eigenvalue (2) Fourth Order Systems 
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Figure 4.8. Unit Impulse Response of the Compensated (1) and 
Prescribed Eigenvalue (2) Fourth Order Systems 
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Figure 4.9. Unit Step Response of the Compensated (1) and 
Prescribed Eigenvalue (2) Fourth Order Systems 
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Figure 4.10. Unit Impulse Response of the Uncompensated (1) 
Prescribed Eigenvalue (2) Seventh Order Systems 




117 



Figure 4.11. Unit Step Response of the Uncompensated (1) and 
Prescribed Eigenvalue (2) Seventh Order Systems 
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Figure 4.12. Unit Impulse Response of the Compensated (1) an 
Prescribed Eigenvalue (2) Seventh Order Systems 
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Figure 4,13. Unit Step Response of the Compensated (1) and 

Prescribed Eigenvalue (2) Seventh Order Systems 
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V. CONCLUSIONS AND DISCUSSION 



A method has been shown for determination of the state 
weighting matrix in order to satisfy a prescribed set of 
eigenvalues through phase variable state feedback. From a 
strictly mathematical viewpoint, this technique requires only 
a knowledge of matrix algebra. Every attempt has been made 
to avoid the necessity of inverting a matrix. The intro- 
duction of Chapter IV made known the fact that previous devel 
opments in this area have suffered the main drawback of res- 
triction. The author believes the method presented here, 
using Cauer Second Form, overcomes many of these restrictions 
It presents a rational computational procedure for determina- 
tion of the weighting matrix, Q; the system eigenvalues are 
only required to be in the left half of the complex plane as 
opposed to the left of a line parallel to the imaginary axis ; 
and the method is no more complicated for multiple or complex 
eigenvalues than a system with linearly independent eigen- 
values (or eigenvectors). 

It should be noted that some authors ’’define" the 
eigenvalue ( s ) of a matrix to be only the real root(s) of the 
characteristic equation. In the development that has pre- 
ceded in this thesis, all roots are considered eigenvalues 
of the associated matrix. 
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The algorithm derived in Chapter IV, is designed as a 
basis for future research. In particular, if the designer 
is working with nth order systems where n is relatively 
large, it may be advantageous to look at using mth order 
simplified models (m<n) of each system by a partitioning 
scheme similar to that in Chapter III. Again, it is em- 
phasized that reduced order models do not necessarily yield 
stable systems. If the simplified system retains the basic 
characteristics of the original system, especially in steady 
state, then this would appear to be a reasonable approach. 

A parallel approach could also be inventigated regarding 
multi-input multi-output (MIMO) systems, treating each 
element of the system transfer matrix as an individual 
transfer function. 

To the author's knowledge, no work has been done in the 
digital or sampled data areas involving continued fraction 
theory. This area should be considered due to the increasing 
use and need for digital control systems . 

These topics represent just a few of the areas available 
for future work. 
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APPENDIX A 



Consider the system 

X = Ax + Bu, (A-1) 

where x is an n-dimensional state vector, u is the input 
function, and A and B are time-invariant (nxn) and (nxl) 
matrices, respectively. The phase variable (canonical) 
system representation is defined as 

/V 

V = Av + Bu, (A- 2) 

W .s, .S, .s, 

where v is an n-dimensional state vector and 





0 1 0 0 




0 




0 0 1 0 




0 




0 0 0 0 




0 


/\ 




/\ 




A = 


• • • • 


, B = 


• 




• • • • 




• 

• 




0 0 0 1 




• 








1^ 




_jl 2 3 nj 







(A-3) 



The systems represented in (A-1) and (A-2) are said to be 
equivalent if and only if there exists a non-singular matrix, 
K, such that 



X = Kv 



(A-4) 
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Kalman [23] has shown that a necessary and sufficient condi- 
tion for such an equivalence to exist is that the system 
in (A-1) be completely controllable. 

The controllability matrix of system (A-1) is defined by 



E = [B 


AB : A^B 


• • • • 


: A^“^B] 


(A-5) 


•s. •>» 


• 









or in an equivalent manner 

E — Ce2^e2*«*.s^]5 (A— 6) 

where the (nxl) vector e^ is recursively defined as 

e^^j^ = A e^ , e^ = B . (A-7) 

The controllability matrix of system (A-2), E, is defined 

/N ^ 

in a similar manner with A and B. Since there is only one 

control input, a necessary and sufficient condition for 

/\ 

controllability is that the (nxn) matrix E (or E) have an 
inverse . 

Silverman [20] has shown that if the system in (A-1) 
is controllable, then the transformation matrix, K, is 
determined by 

K = E E~^ , (A-8) 

where 
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<1 



{ 






n 



1 

0 



S ^ 

^5 

^6 

0 

0 



n 

0 

0 



0 

0 



(A-9) 



and 



a = 



n-1 

a 

n 



= -E"^ B, 



(A-10) 



The elements of a are the coefficients of the characteristic 
polynomial : 



/N ri • -1 

detCSI-A] = detCSI-A] = S + I a.S^"^ 

^^1 1 



(A-11) 



The matrix inversion in equation (A-10) can be avoided by 

using the Leverriexi-Fadeev method for calculating the 

coefficients of the characteristic polynomial. Once the 

•^_1 

coefficients are known, E is written by inspection. 
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Rane [24] presented a simplified procedure for finding 
the transformation matrix, K, requiring no matrix inversions. 
Substituting equation (A-4) into (A-1) and premultiplying 
equation (A-2) results in: 



X = AKv + Bu 



= KAv + KBu 



(A-12) 

(A-13) 



Comparison of equations (A-12) and (A-13) yields 



AK = KA 



(A-14) 



and 



B = KB . 



(A-15) 



Partition K into n column vectors, each (nxl), so that 



K = [k, .... k ] . 



(A-16) 



Substitution of (A-3) and (A-16) into (A-14) and (A-15) gives 



A[k, k„ k, . . . . k ] = 
^ i ^6 ~n 



[k, k^ k„ k ] 

~ 1 ~ z ~6 ~n 



0 

0 



1 0 0 

0 0 0 



-a^ -a2 ^n-1 “^n 



(A-17) 
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and 



B = [Ic, k- . . . . k ] 
~ -vl ~2 ~3 ~n 



0 

0 




(A-18) 



1 



From (A-17) and (A-18) 



k = B 
~n 



k„-,=Ak + k a 

~n-l ~ ~n ~n n 



k o=Ak T+k a , 
~n-2 ~ ~n-l ~n n-1 



k« = A k^ + k^ a- 

^ i, -w-wo '^n 0 



Si = ^ !52 * !$n ^2 ’ 



(A-19) 



or, in general, 



k ,=Ak -.T+k a 
~n-i ~ -n-x+1 ~n n-i+1 



(A-20) 



for icCl, 2, ..., n-1]. The column vectors k^, ..., k^ are 
found in a simple recursive manner and completely 'determine 
the transformation matrix. 
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APPENDIX B 



INVERSION OF CAUER I AND CAUER II FORMS 



This program was written in FORTRAN IV, and requires 
minimal input . The only information required is : 

1. the order of the system 

2. which inversion is required 

3. the quotients from the continued fraction expansion. 
Multiple data sets are possible, and input in the following 
format : 



Card 


Columns 


Description 


Format 


1 


1-3 


M = the desired order transfer 
function 


13 


2 


1-20 


h, , the first quotient of either 
Cauer I ;or.. Cauer II continued 
fraction expansion 


D20.13 




21-40 


h 25 the second quotient 


D20.13 




41-60 


hg , the third quotient 


D20.13 




61-80 


h|, , the fourth quotient 


D20.13 


• 

• 

t 


• 

• 

• 




• 

• 


N 


1-20 


7 ’ (4N-7)th quotient 

sdch^that 4N-7<2M 


D20.13 




21-40 

41-60 

61-80 


Nn-6’ ‘+N-6<2M 
Nn-5’ ‘+N-5<2M 

h4N-4’ 4N-4<2M 


D20.13 

D20.13 

D20.13 


Four quotients per card until 
2*M quotients have been input, 
where M is the system opder. 
Assume this is the Lth card. 
The (L+l)th data card begins 
the second data set. 
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L+1 



1-3 

4-6 



13 

13 



M=system order 

K=1 for Cauer I and K=2 for 
Cauer II inversion 



The computer program has been written to handle up to 
20th order systems (M_<20). If it is required to work with 
higher order systems, only one card change must be made. 

The specification statement is modified to read: 

REAL*8 A(N), B(N), C(N)/N*0.Z, D(N)/N*0./, DZERO 

where N is an integer no larger than 999. This restriction 
can be lifted by changing statement 2 to read: 

2 F0RMAT(2IR) 

where R is the mantissa of log^^CN). The REAL*8 in the 
specification statement indicates that all following varia- 
bles and arrays are real valued and in double precision. 
Modification to either single or extended precision would 
require changes in all format statements. If this is 
desireable, the user should consult references [25] and [26].- 
Execution time has shown to be less than .18 seconds 
for systems of order 10 or less. 
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// EXeC FORTCLG 
//FORT.SYSIN DD » 

REAL*8 A(20) ,B(20) , C( 20)/20*0 ./ , D (20) /20*0 ./ , DZERO 
... READ IN SYSTEM ORDER, a WHICH CAUER INVERSION ... 

16 REA0(5,1) M,K 

1 F0RMATC2I3) 

... READ IN QUOTIENTS FROM C3NT. FRAC. EXPANSION ... 

REA0(5,2) (A( I) ,B( I ) ,I=l,M) 

2 FORMAT! AD20. 13) 

... DETERMINE IF INVERSION IS CAUER I OR CAUER II ... 

IF (K.EQ.l ) GO T3 10 

... CAUER II INVERSION ... 

... INITIALIZATION ... 

C(M) = B(M) 

Ml = M-1 
D(M1) = 1.0 
D(M ) = A(M)*C(M) 

... ITERATION ... 

DO 3 1=1, Ml 
L = I + l 
K = N-L 
DO 4 J=1,L 
KJ = K+J 
KL = KJ-1 

[CO = I 

C(KJ) = BiKI^DUJ) 

IF(J.NE.L) C(KJ) = C(KJ)+C(KP) 

IF(KL.EQ.O) GO TO 5 
0(KL) = A(K)*C(<L) +D(KJ ) 

GO TO 4 
5 DZERO = 1.0 
4 CONTINUE 

DIM) = A(K-H)*C(M) 

3 CONTINUE 
GO TO 11 

... CAUER II INVERSION ... 

... INITIALIZATION ... 

10 C(M) = B(M) 

MM = M-1 

D(MM) = A(M)*C(M) 

D(M) = 1.0 

... ITERATION ... 

DO 6 1=1, MM 
IP = I+l 
MI = M-IP 
00 7 J = l, IP 
MJ = MI-^J 
ML = MJ-1 
MP = MJ+1 

C(MJ) = B (MI ) #D(MJ )+C(MJ) 

IFIML.EQ. 0) GO TO 8 
0(ML) = AIMD^CIMJI+DIML) 

GO TO 7 
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8 OZERO = 1.0 
7 CONTINUE 
6 CONTINUE 

... WRITE OUT TRANSFER FUNCTION IN TERMS OF 
NUMERATOR AND DENOMINATOR COEFFICIENTS 
WITH APPROPRIATE POWER OF S ... 

11 WRITE(6,12) 

12 F0RMAT(///11X, 'NUMERATOR* , I 5X , ' DENOM INATOR • , 
llOXt 'POWER OF S* ) 

WRITE(6,13) OZEROfM 

13 FORMAT(///31X,D20.13tlOXt 12) 

00 14 I=l,iM 

MM I = M— I 

WRITE(6,15) C(I),D( I),MMI 

14 CONTINUE 

15 F0RMAT(///6X,02D.13,4X, 020 . 13 , lOX , 1 2 ) 

GO TO 16 



FOR ACTUAL RUN THIS CARO IS /♦ IN COLUMNS 1 ANO 2 
.SYSIN 00 * 

OATA INPUT 
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APPENDIX C 



DETERMINATION OF WEIGHTING MATRIX 
(Q) FOR PRESCRIBED EIGENVALUES 

This FORTRAN IV program was used exclusively on the 
Naval Postgraduate School’s IBM 360/67 digital computer 
and includes the associated job control language statements. 
The program consists of a main program and nine subroutine 



subprograms . 
below. 


The purpose of each subroutine is delineated 


SUBROUTINE 


DESCRIPTION 


READ 


read in coefficients of numerator and denomi- 
nator polynomials of both transfer functions , 
and places each system in phase variable form. 


RAMAT 


determines the Routh array matrix and the 
transformation matrix, P. 


MULT PH 


multiples two matrices, Y and Z, and gives 
the resulting matrix YZ." 


HMATRX 


determines the quotients of continued fraction 
expansion, H1(H2), and the state matrix in 
Cauer II fofm,''HHl(HH2) . 


POLYNM 


T 

determines the product det(SI-A)xdet(SI+A.. ) 


HELP 


computes the matrices G, T, and W as given 
in Chapter IV. 


QTILDA 


computes the matrices, Q, and Q from results 
of subroutine HELP. 


QPIJ 


determines the^off diagonal elements, q^^ > 
of the matrix Q. ^ 


WRITE 


writes all two-dimensional matrices. 
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The input required has been reduced to a minimum. Multiple 
data sets are possible; and input as indicated below: 



Card Columns 



Description 



Format 



2 1-16, 
17-32, 
33-48 



Lf2 1-16, 
17-32, 
33-48 



2L+2 1-16, 

17-32, 
33-48 , 



3L+2 1-16, 

17-32, 
33-48 , 



N; the system order 

a for ieCl , 2 , . . . , n] ; 
tRe'^denominator coefficients of 
the known transfer function, 
for an Nth order system, L 
cards are required, where L 
=K+1 and K is the integral part 
of N/5. 

b^_^ for leCl, 2 , . . . ,n] ; the 
numerator coefficients of the 
known transfer function, L 
cards required. 

a for iiCl,2, . . . ,n] ; 
tRe'^denominator coefficients 
of the transfer function with 
prescribed eigenvalues. 

3 , for isCl, 2 , . . . ,n] ; the 

numerator coefficients of the 
transfer function with pre- 
scribed eigenvalues. 

for multiple data sets, repeat 
the same prodecure. Each data 
set requires 4L+1 cards. 



13 

F16 .8 



F16 .8 



F16.8 



F16.8 



This program has been written to accept systems up through 
20th order. To increase the capability of the program, only 
the dimension statements and the second continuation card of 
the equivalence statement require modification. The system 
order capability can be increased to 50. Beyond 50th order, 
the program requires an excessive amount of storage space 
(>510K bytes). Even this limitation is easily overcome by 
removing the four cards between statements 1000 and 1001. 
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P 



Any other modifications (i.e., single or extended pre- 
cision) require extensive changes to all subprograms. In 
this case, the user should consult [25] and/or [26]. It is 
recommended that an object deck or disk storage be used when 
available as compilation time is approximately 70-80% of 
total CPU time. 
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SUBROUTINE POLVNMIN, DB, E» 

IMPLICIT REALMS (A-H,0-Z», INTEGER (I-Nf$ 
DIMENSION 0B(20I, E(21lf 0Cl(21)t 0C2(2U 



fM 



O 

Q 

* 



O 

Q 


















^0- 








CO 
















f-4 














o 


r- 




CM 


o 


ou. 






o 








o 


h— —J 






•f CO 


o 




CO 


+ 


< 






Q •• 


#^i^ 




►- 


Q 


ox 












m4 


0^-4 








• 1— 


wO 






•l-l- 


•w* 










o 


l-Z 


ozz 


-i^LU 






••1— lO 


Q 


Z ••CO 


•» •» 








II Q-i- 


1 


*rsJN. II ^-0 


• II 




^ II 


II -5 


> m4 




CO II 


II II ZLU 




+ II UJ 


1 


II « 


II UJ 




—-iz+^z— ujoiai 




2 O 




Gf^ 


1 rD 


II II u. 


-j •u.r>r33 




0--2 




-^lU 1 


-^z 




-liHvO 


X-IZZZ 


z 


II 




C0^1L<^^ II 


a: 


r-wh- 


wfN- 




^1- 






o 




^ II CO— II cgz 


< 




X«^*-»ZZZ 


KQ 


H-OOO 


oo 


OUL 


oo 


l-OX 


WOOCLU.WOOO 


UJZ 


ZQQO 


OO^t^Oi-iUOO 


za>»uJoa^*-iujooo 


C^Ui 


O 






m4 


CO 




v0lTl>t' 












•— 1 




^0^4 ^4 





















O OOO O O O 



141 



SUBROUTINE HELP(N) 

IMPLICIT REAL*8 (A-H,0-ZI, INTEGER (I-Nt$) 
DIMENSION H(20,20), AID(20,20), V(20,2D), F(20,20, 
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MIPl = N-IEYE<-1 

P(IEYE,IEYE) = (-C0UNT«-E2( IEYE)-El(IEYEn/(CD(IEYE,NMIPl)* 



IF(IEYE,NMIPI,IEYE|} 
CALL QPIJ(NflEYE) 
311 CONTINUE 
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